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We study a two filament driven lattice gas model with oppositely directed species of particles 
moving on two parallel filaments with filament switching processes and particle inflow and outflow 
at filament ends. The filament switching process is correlated such that particles switch filaments 
with finite probability only when oppositely directed particles meet on the same filament. This 
model mimics some of the coarse grained features observed in context of microtubule (MT) based 
intracellular transport, wherein cellular cargo loaded and off-loaded at filament ends are transported 
on multiple parallel microtubule (MT) filaments and can switch between the parallel microtubule 
filaments. We focus on a regime where the filaments are weakly coupled, such that filament switching 
rates scale inversely as the length of the filament. We find that the interplay (off)loading processes 
at the boundaries and the filament switching process leads to some distinctive features of the system. 

These features includes occurrence of variety of phases in the system with inhomogeneous density 
profiles including localized density shocks, density difference across the filaments and bidirectional 
current flows in the system. We analyze the system by developing a mean field (MF) theory and 
comparing the results obtained from the MF theory with the Monte Carlo (MC) simulations of the 
dynamics of the system. We find that the steady state density and current profiles of particles and 
the phase diagram obtained within the MF picture matches quite well with MC simulation results. 

These findings maybe useful for studying multi-filament intracellular transport. 


I. INTRODUCTION 


One-dimensional driven diffusive systems, unlike their 
equilibrium counterparts are known to exhibit boundary 
induced phase transitions (THU- Such systems have also 
served the purpose of providing a framework for studying 
wide class of driven biological phenomenon ranging from 
transport across biomembranes [4], to transport on in¬ 
dividual cellular filament and cytoskeletal filament 
network El- 

Filament based intracellular transport involves oppo¬ 
sitely directed motors, which use multiple arrays of cy¬ 
toskeletal filaments, to actively transport cellular cargoes 
such as mitochondria, endosomes, and pigment gran¬ 
ules 0 S3- It has been observed that long-distance 
cellular cargo transport on microtubule (MT) filaments 
is achieved by sets of oppositely directed motor proteins, 
e.g; dynein and kinesin , which attach to the cellular car¬ 
goes and transport them actively along these filaments. 
The transport itself is determined by different processes 
at play at the molecular level, e.g; the motor processiv- 
ity, directional switching dynamics of the cargo carried 
by the motors, the underlying filament organization, the 
(un)binding characteristics of the motors to the filament 
and the boundary input (output) rate of cargoes at fil¬ 
ament ends [13] . It has also been observed that both 
long distance regulated transport [l3] and phenomenon 


of jamming arise out of the collective action of the these 
motor proteins 0 El- 

One of the approaches to study transport in such sys¬ 
tems has been to describe it in terms of coarse-grained 
driven lattice gas models wherein the MT filament is 
considered as a 1-d lattice, the interactions between the 
transported cargoes are included via excluded volume ef¬ 
fect and the underlying driven stochastic dynamics due 
to various processes incorporated in the description [§}. 
Some of the previous theoretical attempts have focused 
on the interplay of stochastic directional switching mech¬ 
anisms, directional hopping of individual cargoes and 
the effect of the input and output of the cargoes on the 
boundaries on a single filament SIE3- However for 
a variety of biological situations such as axonal trans¬ 
port in neurons, cargo transport takes place on a multi¬ 
ple parallel array of MT filaments [15], [ljj ■ For example, 
in vitro studies on cultured neurons have revealed that 
cargo switching between neighbouring filaments occurs 
in axonal transport of mitochondria on neurons [18] . On 
theoretical grounds, it has been argued that even with¬ 
out considering the effects of the boundaries, the inter¬ 
play of the translation process on filaments and the fila¬ 
ment switching processes can manifest in form of a phase 
transition between a an inhomogeneous jammed phase of 
the transported cargoes and a freely flowing phase with 
homogeneous cargo density in each filament [l9[. Thus 
studying the role of multiple filaments in determining the 
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transport properties of such systems is of considerable 
importance. 

Driven transport on parallel lattices have been stud¬ 
ied theoretically in different contexts jTol - f^ . and the 
particle switching dynamics between adjacent lanes have 
also been taken into account explicitly in some cases [H- 
[ 22 I . [271 - I30I . [32} . In this paper we will focus on how the 
transport along two parallel filaments is affected by the 
interplay of boundary inflow and outflow of particles at 
the filament ends and filament switching dynamics of par¬ 
ticles. Before we proceed describing the mode in detail, 
we wish to highlight a few aspects of transport that have 
been observed in the context of intracellular transport 
: (a) Experimental studies, such as the one on endo- 
somal transport on MT reveal that cellular cargoes can 
switch between neighbouring filaments [HI, HU. (b) Ex¬ 
periments suggest that cellular cargoes traveling in oppo¬ 
site direction on the same MT can also cross each other 
and continue with their translational motion along the 
same filament (33| . (c) Cargo transport is also depen¬ 
dent on loading and offloading of cellular cargoes at the 
filament ends [lSJ, [34| . Motivated with these experimen¬ 
tal observations, in the model that we study, we focus on 
the interplay of the boundary driven processes of particle 
input and output at filament ends with the active motor 
driven cargo translocation on the filaments and filament 
switching processes. Accordingly, we will consider two 
parallel filaments with oppositely directed species of par¬ 
ticles which translate on the lattices with a specified hop¬ 
ping rate. The oppositely directed particles are also al¬ 
lowed to pass through each other with a certain specified 
rate on the same lattice. The particles are also allowed to 
switch between the lattices with certain finite probability 
only when oppositely directed species meet each other on 
the same lattice, so that the switching between the lane 
is a correlated process [l9|, HU . Thus implicitly we will 
consider that individual motors carrying the cargo have 
a propensity to switch between different filaments when 
they experience a force, when hindered by an oppositely 
directed particle moving on the same filament. This in 
turn can decrease the binding affinity of the motor to the 
filament and induce it to switch and bind to the neigh¬ 
bouring filament. Finally we will allow particles to enter 
and exit the filament ends with prescribed rates. 


density profiles for the particles and construct the MF 
phase diagram for the system and compare these results 
with the Monte Carlo (MC) simulation results. Finally 
in Section [V] we discuss these results in context of multi- 
filament intracellular transport. 


II. THE MODEL 
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FIG. 1: Schematic representation of the processes of translo¬ 
cation, switching between the two filaments, and entry/exit 
at filament boundaries. 


We consider two parallel cellular filaments, represented 
by two finite and parallel one dimensional lattices of 
length L with N sites, with lattice spacing e = L/N. 
The cellular cargoes transported along these filaments 
will be referred to as particles, and will be characterized 
by two different species. Specifically, along each filament, 
the transported cargo can either be a (+) particle which 
moves from left to right on the filament or a (—) particle 
that moves from right to left. Without loss of generality, 
the two filaments are labeled as 1 and 2. The instan¬ 
taneous state of the system is described in terms of the 
occupation numbers, which indicate the spatial localiza¬ 
tion of the two species of particles on the two parallel 
filaments. Specifically, corresponds to a occupation 
number of a particle at site i moving to the right on fila¬ 
ment 1. The maximum allowed occupancy at any lattice 
site is 1 so that each lattice site is occupied either by a 
(+) particle, (—) particle or is vacant, (0). 


In Section[TT]we specify the model, the dynamical rules 
on the lattice and the corresponding equations of motion 
for the system. In Section IIIII we set up the Mean Field 
(MF) continuum equations for the system and analyze 
the boundary conditions at the ends of the lattice. In 
Section m we analyze the different phases possible for 
the system, obtain the corresponding MF steady states 


The dynamics of this system can be expressed in terms 
of the movements allowed for the particles. For the sites 
in the bulk , in each individual filament, a (+) particle can 
hop from site i to site i + 1 with a rate a if that site is 
vacant. Similarly a (—) particle can hop from site i to 
site i — 1 with an identical rate a if it is vacant. If the 
i-th site on filament 1 is occupied by a (+) and if the 
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neighboring site to the right, i.e; site i + 1 is occupied by 
an oppositely directed (—), then two different processes 
can occur: with rate /3 the two particles can swap their 
positions, and the (+) particle moves to the site i +1 and 
the oppositely moving (—) particle moves to site i; while 
with rate tt° 2 the (+) particle can switch the to filament- 
2, at the same corresponding site, identified by its index 
*, if the site is vacant. Similarly, a (—) particle from 
filament 1 at site * +1 can switch to the site of filament 2 
with the same site index i+ 1 with a rate /i° 2 if that site on 
the other filament is vacant. Identical processes that we 
have described for filament 1 also happens for filament 2. 
The rates of filament switching processes from filament 2 
to filament 1 for the (+) and (—) particles are 7r 21 and p 21 
for (+) and (—) particles, respectively. Particle switching 
between the two filaments can be understood as arising 
from the stronger loading force experienced by the motor 
proteins which carry the cargo when they push against 
an oppositely directed cargo. This leads to an increase in 
the rate of motor detachment (along with the cargo) from 
the filament and offers the possibility of a subsequent 
reattachment of the molecular motor to the neighbouring 
site on the other filament [35]. Since we are interested 
in the regime where the filaments are weakly coupled, we 
probe the regime where the filament switching processes 
compete with loading and offloading processes at filament 
boundaries. We systematically implement it by choosing 
the filament switching rates at individual lattice site such 
that they scale inversely with system size so that we have, 

—O _ 7Q2 ,,0 _ 12 —O _ IT 21 ny. J ,.0 _ U21 

^12 — N > "12 — N ’ 7r 21 — N clIiU ^21 — N ’ 

Although the bulk processes are analogous to those in¬ 
troduced in Ref. |19|, [20| , which focused on the collective 
behavior of cargoes moving along filaments with a closed 
ring morphology and with overall particle number con¬ 
servation, we will concentrate on the behavior of such 
particles for open filaments. In this configuration the 
overall particle number is not conserved, and the motion 
on incoming and outgoing particles at the filament ends 
must be accounted for. 

For the boundary sites at the filament ends, a (+) par¬ 
ticle can enter the left end of the filament (with site label 
* = 0) of filament 1 with a rate if it is vacant, and 
it can leave from the right end of the filament (with site 
label IV — 1) of filament 1 with a rate . Similarly a (—) 
particle can enter the right end of the filament 1 with a 
rate if it is vacant and it can leave from the left end of 
the filament with a rate /3-f. Similar processes also occur 
in filament 2 with the corresponding rates being a 2 , P 2 , 
a 2 and jd 2 respectively. All the the dynamic processes 
that characterize the model are schematically depicted in 
Fig. 1. 


III. MEAN FIELD EVOLUTION EQUATIONS 


The time evolution for the average occupation number 
for the two oppositely directed species along each indi¬ 
vidual filament can be expressed in terms of gain and 
loss terms arising from translation and filament switching 
processes. As described in Ref. [2C), these terms involve 
the averages of the local occupation numbers for the par¬ 
ticles at each site as well as various combination of aver¬ 
ages of two-point correlators of site occupation numbers, 
which account for the role of particle correlations in the 
particle collective dynamics. 

We set a = /3 = 1 and choose 7 Ti 2 = 7T2i = H 12 = p-21 = 
7r, which correspond to a symmetric scenario where the 
propensity to switch filaments is the same for both (+) 
and (—) species and it is symmetric about filament 1 and 
filament 2. The Mean Field (MF) evolution equations 
are obtained factorizing the two point correlators of the 
occupation numbers. The continuum MF evolution equa¬ 
tions are derived by rescaling the total length L to 1 and 

letting N —> 00 so that e —>■ 0 @,E3- 

Correspondingly, pi{x), p 2 (a;), n\(x) and ?z 2 (x) are 
then the average densities as a function of the relative 
position in the filament, x. The MF continuum equa¬ 
tions in the bulk can be expressed as, 


d t pi = en [p 2 n 2 (l - Pi - ni) - pini(l - p 2 - n 2 )\ 

- ed x [pi(l-pi)] + 0(e 2 ) (1) 

d t p 2 = £tt [pini(l -p 2 - n 2 ) -P 2 n 2 {l - Pi - «i)] 

- ed x [p 2 (l -P 2 )] + 0(e 2 ) (2) 

d t ni = 67T [p 2 n 2 (l -pi - ni) - pim(l - p 2 - n 2 )\ 

+ ed x [^i(l — ^ 1 )] + 0(e 2 ) (3) 

d t n 2 = en \pini(l - p 2 - n 2 ) - P 2 n 2 (l - Pi - ni)] 

+ ed x [n 2 (l - n 2 )] + 0(e 2 ) (4) 

where we have displayed terms up to first order in e. The 
corresponding expression for the currents of each species 
read, 

J i = Pi(l-Pi) (5) 

^ = -ni(l-ni) (6) 

J 2 + = P 2 O--P 2 ) (7) 

J 2 = —n 2 (l — n 2 ) (8) 
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A. Steady State profiles 


From Eqs. ©-© , the steady state profiles correspond¬ 
ing to the continuum MF evolution of the molecular mo¬ 
tors satisfy, 


dJ+ 

dx 

dJi 

dx 

dJ 2 + 

dx 

dJ‘2 

dx 


= 7 t[p 2 ^ 2(1 - Pi- ni) P\n\ (1 -p 2 - n 2 )\ 

= Tr[p 2 n 2 (l - Pi - «i) - pini(l - p 2 - n 2 )\ 

= 7 r[pini(l -p 2 - n 2 ) - p 2 n 2 ( 1 - pi - ni)] 

= 7r[pini(l -p 2 - n 2 ) -P2n 2 (l - pi- Til)], 


(9) 

( 10 ) 

( 11 ) 

( 12 ) 


which govern the bulk profiles of (+) and (—) particles in 
the two filaments. One can rewrite them in terms of the 
fluxes of the total number of particles and its difference 
introducing 


Jp =Pi(l -pi) + P 2 O- -P 2 ) (13) 

Jn = ni(l - m) + n 2 (l - n 2 ) (14) 


Since in the bulk the total current for the combined sys¬ 
tem comprising of the two filaments is isolated, J p and 
J n have to be spatially constant. Subtracting Eq. m 
from Eq. © we identify an additional conservation law, 

Ji=pi(l-pi) + ni(l-ni) (15) 


which corresponds to the sum of the absolute values of 
the currents of opposite species along one filament. As 
shown in Appendix B, it is useful to reorganize the three 
independent conserved quantities, J p , J n and J\ and 
which remain spatially uniform in the bulk, in terms of 
three new parameters, J 2 = J p — J\ + J n , C\ = ,J P — J 2 
and C 2 = J p — J\. Using these new conserved quanti¬ 
ties, the equations for the density profiles in the bulk can 
be decoupled, and express them in terms of one single 
density variable. Specifically, we can write 

dPl 1 r ± ± ±\ 

HI = n) 

- l4i v Pi (l -Pi -dp,)] (16) 

where 77 ^, and are functions of p\ alone. Their 
explicit functional dependence is provided in Appendix 
B. Similarly, we can get decoupled differential equations 
for m, p 2 and n 2 , 

dfl 1 1 r ±1-1 ± ± \ 

HI = TVlsd’™-. t 1 

- Mn^ni (l-m-r?^)] (17) 


dp 2 

dx 


~ ! J 2 p 2 [ P2iq P » “ Vp 2 - "£) 

^P, U P2 (1-P2-^)] 


(18) 


dn 2 

dx 


1 -2n 2 [ni7? " 2 ^ “ & ~ 

(l -« 2 -»?*)] (19) 


The explicit form of the decoupled differential equations 
which govern the density profiles and the relevant so¬ 
lutions are discussed in Appendix B. In Appendix B, 
Eo. (IBHIB4l) . provides the mathematical expression for 
the quantities present in Eq. mm- Appendix B also 
describes the different sets of density profiles that can be 
obtained from the previous set of equations. We have 
found that there are 16 different branch solutions cor¬ 
responding to the decoupled differential equation. The 
choice of a particular boundary condition corresponding 
to a particular phase, restricts the possible choices to 8 . 
Finally as discussed in Appendix B, physical considera¬ 
tions such as bounds on the physical value of density se¬ 
lects an unique solution to these differential equations for 
each set of prescribed boundary conditions. Therefore, 
the continuous MF equations determine the density pro¬ 
files of the particles in the two filaments, once the bound¬ 
ary conditions are prescribed. Although Eas. (fTUl - [I!)l) for 
the different species densities decouple in the bulk, they 
are coupled through the boundary condition of Eqs. CES 
m- In the next subsection we discuss the set of boundary 
conditions satisfied by these differential equation corre¬ 
sponding to a particular phase. 


B. Boundary conditions 


The allowed phases that characterize the state of trans¬ 
port on the two filaments is controlled by the particle 
input and output at the boundaries. The steady state 
density profiles are determined by first order differential 
equations. This is due to the fact that the diffusive con¬ 
tribution is of higher order in the lattice spacing, e, and 
their contribution drops in the continuum limit, e —> 0 . 
As a result, the system cannot fulfill, generically, the in¬ 
put and output boundary conditions and boundary layers 
are expected |<3]. One then must determine which of the 
fluxes at the filaments’ ends determine the bulk density 
profiles and under which conditions coexistence between 
density phases can develop along the filaments. 

For the system composed of two filaments, we have 8 
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different particle entrance or exit rates at the boundaries, 
that we express as ajy 2 ,/3jy 2 . However there are only 4 
boundary conditions to be specified for the complete so¬ 
lution of the steady state differential equations. In order 
to figure out the possible physically relevant boundary 
conditions, it is useful to build upon the boundary con¬ 
ditions that are satisfied by a closely related 1-D lattice 
gas model [I, 3, which has the similar translocation dy¬ 
namics along the filament as our model, but which does 
not allow for inter-filament exchange processes. In fact 
the model discussed in Ref. 1,2] is exactly the same as 
ours for the particular case of ir = 0, which corresponds 
to a situation where the inter-filament switching dynam¬ 
ics of the particles is turned off. Since for our case, the 
filaments are weakly coupled thus it is expected that the 
boundary conditions satisfied for our two-filament sys¬ 
tem are the various possible combination of the bound¬ 
ary conditions that are satisfied for individual lattice for 
the case studied in Ref. 00- However we would like 
to stress that although the boundary conditions are ob¬ 
tained as simple combination of boundary conditions of 
the individual lattices, the resultant density and current 
profile obtained by spatially integrating the steady state 
differential equation would be qualitatively different due 
to the lattice switching term in the bulk. 

We enumerate the possible combination of the bound¬ 
ary conditions and the resultant phases for each of those 
particular combinations. 

Filament 1 in HL phase and filament 2 in HL phase 
((HL) i — (HL) 2 ): When the bulk current of (+) matches 
with the output current of (+) at the right boundary and 
the input current of (—) matches with the bulk current of 
(—) at the right boundary both for filament 1 and 2, the 
resultant phase corresponds to a situation where the (+) 
particles are in high density) H) phase and (—) particles 
are in low density(L) phase (refered to as HL phase), in 
both the filaments. The boundary conditions that are 
satisfied in the continuum limit are, 

J ir = P 1 P 1 R = Pir{ 1 ~Pir) 

J\r = - 04(1 -PiR - Rir) = - nin) 

J 2R = Pi PiR = Pir( 1 ~ Pir) 

J 2 R = ~ a l (1 ~ P2R — R2r) = — n2R(l — U2R.) (20) 

where J+ R , Jfi R , J^ R and J^ R refer to the currents for 
(+) and (—) particles in filaments 1 and 2 at the right(R) 
boundary respectively, while Pir 1 nm, p 2 R and n 2 R refer 
to the densities of (+) and (—) particles in filaments 1 
and 2 at the right boundary. Using Ea. (l20l) . we can find 
the expression of the boundary densities at the right end 


of both filaments in terms of the entry and exit particle 
rate 


PiR = 

niR = 

P2R = 

l-/3i + 

(1 + a i ) ~ 1 

\J{ 1 + aq ) 2 

- 4 a^Pi 

1 -Pi 

2 


n 2 R = 

(1 + a 2 ) — 1 

\] ( 1 + a 2 ) 2 

— 4a 2 P 2 


2 



By symmetry there can be another phase where both in 
filament 1 and 2, (—) are in high density phase while the 
(+) are in low density phase and the current conditions 
are satisfied at the left boundary. Further one can find 
a situation where for filament 1 the current at the left 
boundary for (+) and (—) matches with the bulk current, 
while for filament 2, the current at the right boundary 
for (+) and (-) matches with the bulk current. Similarly, 
there exists a phase for where for filament 1, the cur¬ 
rent at the right boundary for (+) and (—) matches with 
the bulk current, while for filament 2, the current at the 
right boundary for (+) and (—) matches with the bulk 
current. For all these 4 different phases, the structure of 
the boundary condition is exactly similar. 

Filament 1 in LL phase and filament 2 in LL phase : 
((LL) 1 — (LL) 2 ): When the bulk current of (+) particle 
matches with the input current of (+) at the left bound¬ 
ary and input current of (—) matches with the bulk cur¬ 
rent of (—) at the right boundary both for filament 1 and 
2, the resultant phases corresponds to LL phase in both 
the filaments. The boundary conditions satisfied by the 
currents are, 

Jil = af (1 - pil - ni L ) = pi L (l - pil) 

J 1 R = -«r(l -PiR -Rir) = -nifl) 

J 2 L = a 2~(l ~P 2 L - n 2L ) =P2i(l ~P2l) 

J 2 r = — oi 2 (1 — P 2 R — n 2 n) = — n 2 n(l — n 2 n) (22) 

while the expression for the currents at the other bound¬ 
aries read 

J t R = Pt PiR 
= Pi n lL 

J 2 R = Pi P2R 

J 2 L = P 2 n 2L (23) 

Using Eq. (E!?l) and Eq. (P?3l) , we can express the boundary 
densities as a function of entry and exit rates and the 
currents [ 2 ], 
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at (! -P2L -n 2L ) = _ Jzl _ 

1 — P2L J2L/ a 2 + ^2 hi P 2 

P 2 P 2 R _ J 2 R 

l~P2R~ I-J^r/PZ 

a^" (1 — P 2 R — U2r) _ Jf R 

1 — Vl2 R ^2 r! ^2 + ^2r/ fit 

P 2 n< 2L _ J 2 L (24) 

1 -n*L 1 -J 2 JP 2 

As opposed to the continuity of the overall particle fluxes 
at filament’s ends due to particle conservation, J p l — J p r 
and J u l = JyiR , the current on the left and right end 
of an individual filament track will in general differ due 
to particle filament switching. In order to determine the 
densities for the two types of particles in this phase, we 
will assume that the (+) current in left boundary of fil¬ 
ament 1 equals the (+) current in right boundary of the 
same filament, (J^l) = {^ir}- Similarly we use the same 
equality for (—) particles, (Ji R ) = ( Jf L ). Analogously, 
we equate the currents at the left and the right boundary 
in filament 2, J± L = J^r and J^ L = J^r- This is a rea¬ 
sonable assumption because tracks are weakly coupled, 
as has been checked with MC simulation. 


P2L = 

P2R = 

n 2 R = 

n 2 L = 


This fact allows us to obtain an expression for the 
boundary densities in this phase. Explicitly, from Ea. E^l) 
and Ea. (l24l) . the boundary densities for both the lanes 
read , 


1 


1 


Pil = 1 -- Pil) - —n 1R (l - ni R ) 


1 


1 


n\R = 1 - —tPil(1 — Pil) - -n\R{l - uir) 


a 1 


P 2 L = 1-TtP2l(1 - P 2 L ) - - n 2 R) 

<*2 P2 

ri 2 R = 1 - -rrP2L(l -P2l) - -n 2 R(l - u 2 r) (25) 

P2 a 2 

These coupled algebraic equations for each filament track 
can be numerically solved to get the corresponding 
boundary densities from which the density profiles can 
be numerically derived. 


Filament 1 in LL phase and filament 2 in HL phase : 
((LL)i — (HL) 2 ): When for filament 1, the bulk current 
of (+) matches with the input current of (+) at the left 
boundary and input current of (—) matches with the bulk 
current of (—) at the right boundary and the bulk current 
of (+) matches the output current of (+) at the right 


boundary and the input current of (—) matches with the 
bulk current of (—) at the right boundary in filament 
2, the resultant phase for the system corresponds to LL 
phase in filament 1 and HL phase in filament 2. 

In this case we again assume the continuity of the fluxes 
separately for the two particle types along filament 1, 
JZl = J? R and Jf L = Jf R . Accordingly, the equations 
for the boundary densities are similar in form to those 
expressed in Eq. (E5l) and the boundary densities for fila¬ 
ment 1 can be obtained numerically as discussed earlier. 
For filament 2, the densities are determined by Ea. (l2lT) . 
Due to the symmetry in the swapping rates between the 
two filament, a second phase with analogous structure is 
feasible,where filament 1 in HL phase and filament 2 in 
LL. There is still a further symmetry associated with 
the HL phase in any of the two filaments, i.e; if the cur¬ 
rent of (—) and (+) at the left boundary matches with 
the currents in the bulk, the structure of the boundary 
conditions remains unaltered. In order to illustrate this, 
consider HL phase in a particular filament, then the cor¬ 
responding densities are determined by Eq. m and the 
boundary condition at x = 1 is satisfied, so that (+) is in 
high density phase and (—) is in low density phase. But 
analogously we could have a situation where the bound¬ 
ary condition at x = 0 is satisfied with (+) particles in 
Low density phase and (—) in high density phase. This 
situation would correspond to a different overall phase, 
but structure and the form of boundary density would be 
the same as Eq . (I2T1) with the only difference that the in¬ 
dexes of (+) and (—) in the expression for the boundary 
density in Eq , (l2ll) is interchanged. Thus this structure 
of boundary conditions would correspond to 4 distinct 
phases. 


IV. PHASES AND PHASE DIAGRAM 


Since the boundary fluxes control the particle fluxes in 
the bulk, once we have determined the expression for the 
densities at the filament boundaries in terms of the input 
and output rates at the filament boundaries, the MF den¬ 
sity and current profiles in the bulk can be determined 
from Eq. IHGHlUl) . However, unlike the case of model in 
Ref. pj, i where the steady state density profiles and the 
corresponding phases in the bulk are solely determined by 
the boundary fluxes, the density profiles and phases now 
emerge from the interplay of the boundary processes and 
particle filament switching. This distinctly alters the na¬ 
ture of the density profiles and the topology of the phase 
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diagram for the system under study. 

Specifically, a first major consequence of particle ex¬ 
change between filaments is that the density and current 
profiles in the bulk are no longer spatially homogeneous. 
In fact for certain range of input and output particle 
rates, the competition between the bulk and the bound¬ 
ary processes can result in density shocks along the fila¬ 
ments which are localized in the bulk. 

Moreover, particle change between filaments also al¬ 
lows for phase coexistence in the bulk apart from the 
pure phases which satisfies only one set of boundary con¬ 
ditions. This feature of phase coexistence occurs only 
because the current profiles in the bulk for a set of bound¬ 
ary conditions are not homogeneous so that different MF 
solutions for the current, intersect each other at specific 
spatial location in the bulk of the system. The phase 
coexistence in the bulk happens when the currents of 
the different MF solutions arising out of different bound¬ 
ary conditions match at a point in the bulk of the two- 
filament system. In that case part of the system obeys 
one set of the boundary conditions while the other half 
obeys another set of boundary condition and these set 
of solutions are joined by the condition of continuity of 
current at a particular location in the bulk of the lat¬ 
tice. In general the system selects the set of MF steady 
state solution for which the corresponding current is min¬ 
imum. This holds true as long as the MF current profiles 
do not attain the maximal current value in the bulk. In 
this paper we have restricted our analysis to the region 
of parameter space of entry and exit rates of particles for 
which the condition for maximal current is not reached. 

As a result of these new features, the topology of the re¬ 
sulting phase diagram changes qualitatively with respect 
to the collective behavior in the absence of inter-filament 
particle exchange. In the following subsection we first de¬ 
scribe the procedure to find the MF density and current 
profiles and determine the resultant phase. Subsequently 
we discuss the topology of the resultant phase diagram 
for this system and obtain the equations for the phase 
boundaries. 


A. Density profiles and emerging phases 


In order to find the MF density and current profiles 
using Eg. (11 611191) . we have to first determine the three 
independent conserved currents in the system e.g; J p , J n 


and J\ . Subsequently, we have to provide the appropriate 
values of the densities at the boundaries to completely 
specify the solutions for the individual species. 

(HL )i — (HL )2 phase: Here, first we determine the val¬ 
ues of the boundary densities at the right end of both fil¬ 
aments, e.g; pm, nm, P 2 R and n 2 r using Ea. (l2lT) . Thus 
J p , J n and J\ can be determined at the right boundary. 
The entire density profile can now be found out by evolv¬ 
ing the MF solution from the left end of both filaments 
using Eg. (11611191) . In Fig.2, we show the comparison of 
the MF profile with the MC simulations for this phase. 
It illustrates that both the density and the current pro¬ 
files in the two filaments are not spatially homogeneous 
in contrast to similar phases in Ref. 00- A similar 
procedure can be used to find out the profiles for the 
corresponding (LH) 1 — (LH ) 2 phase. 

(LL) 1 — (HL) 2 phase: The values of the boundary den¬ 
sities, pil and niR, can be determined by numerically 
solving Eo. flUSl) for the first filament. For the second fil¬ 
ament we use Eq. m to determine p 2 R and nm- Thus, 
both J n and J 2 are identified at the right boundary. To 
determine J p we use the method of successive iterations. 
In the first iteration we set pm = pil and obtain J p 
to get the density profiles. Again these profiles are not 
uniform due to inter-filament switching processes. Con¬ 
sequently, the density value obtained at x = 0 by evolving 
Eq. (fTBl) from x = 1 is not same as pm- Hence, we take 
the difference between these two values at x = 0, and this 
difference is added to pm and evolve it to get the new 
density profile. We repeat this process until convergence 
is reached [36). This procedure allows then to derive the 
entire density profile by evolving the MF solution from 
the right end of both filaments using Eq. (11 611191) . An 
analogous procedure is used to identify the profiles for 
the case of (LL) 1 — (HL ) 2 phase, (LL) 2 — (LH) 1 phase 
and (LL) 1 — (LH) 2 phases. 

(LL) 1 — (LL ) 2 phase: For this phase coexistence, we 
need to determine the values of the boundary densities 
Pil, riiR, P 2 L and n 2 R by solving Eq. G5l) and use these 
values to extract the associated boundary fluxes, J p and 
J n . To this end, we assume nm = ri\R, a symmetry that 
holds in the absence of particle filament exchange 0. 
This relation allows us to determine J\ and consequently 
fix the value of n 2 L ■ Again, we use the process of succes¬ 
sive iterations for determining both nm and n 2 L- The 
entire profile can be found out by evolving densities from 
the left end of both filaments using Eos. (116111111) (See 
Fig-3). 










FIG. 2: Steady state (a) density (p) and (b) current ( J) profile 
for (+) and (—) species in filament 1 and 2 as function of 
normalized distance (X) when the system is in (HL)i — (HL) 2 
phase. Here, af = 0.8, oft = 0.2, fit = 0.25, ft" = 0.7, 
ft = 0.9, eft = 0.4, ft" = 0.3, ft" = 0.3, and 7ro = 1.0. 
Points are obtained by MC simulations done for system size 
of N = 2000. Solid lines are the corresponding MF solutions. 
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FIG. 4: Steady state (a) density and ( b ) current profile for (+) 
and (—) when the system is in (HL) i — (SL) 2 phase. Here, 
ft = 0.45, aft = 0.2, ft = 0.25, ft" = 0.7, at = 0.2, a J = 
0.2, ft" — 0.2, ft" = 0.6 and n = 1.0.Points are obtained by 
MC simulations done for system size of N = 2000. Solid lines 
are the MF solutions corresponding to (HL) 1 — (HL) 2 and 
(HL) 1 — (LL) 2 phases. 



FIG. 3: Steady state (a) density and ( b ) current profile for 
(+) and ( —) when the system is in (LL) 1 — (Lift phase. Here, 
at = 0.2, eft = 0.4, = 0.5, ft" = 0.5, ft = 0.7, oft = 0.5, 
fit = 0.4, fix = 0.3, and 7r = 1.0. Points are obtained by MC 
simulations done for system size of X = 2000. Solid lines are 
the corresponding MF solutions. 


(HL) 1 — (SL) 2 : In this phase while filament 1 is in 
HL phase, there is phase coexistence in the second fil¬ 
ament, such that at the right end of the filament, the 
boundary condition corresponding to HL phase is sat¬ 
isfied while at the left end of filament 2, the boundary 
condition corresponding to LL phase is satisfied and this 
phase is characterized by a density discontinuity of (+) 


species in filament 2 which results in a shock profile in 
that filament. For phase coexistence, the profile has to 
be such that the current for (+) corresponding to HL so¬ 
lution equals the current for (+) for the LL solution at 
a particular spatial location between the filaments ends. 
This current continuity condition follows from the fact 
that in the bulk the total added current of the two fil¬ 
aments has to be conserved as there is no particle ex¬ 
change with the environment. The phase coexistence 
can be thought of as a mixture of the (HL) 1 — (HL )2 
phase with (HL) 1 — (LL) 2 . Accordingly, as one changes 
the parameters corresponding to input and output rate 
of particles, starting from a pure (HL) 1 — (LL) 2 phase, 
the system can evolve into (HL) 1 — (HL) 2 through an 
intermediate phase coexistence region corresponding to 
the (HL) 1 — (SL) 2 phase in parameter space, where an 
incipient shock of (+) particles originating at the right 
end of the filament eventually reaches the left end of the 
filament on change of parameters in the phase diagram. 
In order to determine the density profiles for this case, 
we make us of the fact that we know pm, ni_R, p 2 R and 
n 2 R because the right end of both filaments are in the 
HL phase. This property allows us to identify the three 
independent conserved currents e.g; ft, ft, and J n for the 
entire filaments and the density profile from the right end 
of the filament can be plotted using Eqs. (I16H19H . For the 
left end of filament, having determined p\L from Eq. 1251 
the remaining densities at the left boundary can deter¬ 
mined using the three conserved currents. Now the LL 
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profile can be simply determined starting from the left 
end of the filament, using Eqs. (116M19I) . The spatial lo¬ 
cation in the bulk for which the current for this solution 
matches with the current for the other solution obtained 
for HL—HL phase defines the position of shock, as shown 
in Fig.4. 

(LL )i — {SV )2 : In this phase while filament 1 is in LL 
phase, there is phase coexistence in the other filament. 
For filament 2, at the right end the boundary condition 
corresponding to HL phase is satisfied while at the left 
end, the boundary condition corresponding to LL phase 
is satisfied. In order to determine the density profile for 
this case, we use Eq . (12ol) to identify n\R at the right end 
of filament 1. Analogously, we also know p2R and n 2 R 
because filament 2 is in the HL phase. At the left end 
of both the filaments are in LL phase so that p\L and 
P2L are known using Eq .(ESI). Therefore it follows that 
J 2 , J p and J n are known for the entire filaments and all 
the densities at both filament boundaries are thus deter¬ 
mined. Now the density profiles originating from both 
left and right end of the filament can be plotted using 
Eqs. (11611191) separately. The point at bulk for which the 
current for both these solution match defines the shock 
position. 
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FIG. 5: Steady state (a) density and (6) current profile for 
(+) and (—) when the system is in ( SL)i — (SL) 2 phase. 
Here, a+ = 0.35, cq = 0.2, /3+ = 0.292, /3~ = 0.7, at = 
0.2, c %2 = 0.2, /31 = 0.2, @2 = 0.6 and n = 1. Points are 
obtained by MC simulations done for system size of N = 5000. 
Solid lines are the corresponding MF solutions for the phases 
(HL) 1 - (HL ) 2 , (HL) 1 - (LL) 2 and (LL) 1 - (LL) 2 . 


gion close to the left ends is a (LL) 1 — (LL ) 2 phase. In 
between a (HL) 1 — (LL )2 phase develops. The two den¬ 
sity shocks in the bulk separate these three regions. For 
the phase region adjoining the right end of the filament, 
the boundary condition corresponding to (HL) 1 — (HL) 2 
phase is satisfied and the boundary densities are deter¬ 
mined by Eq. m and the MF density profiles are ob¬ 
tained by evolving the MF solution from the right end 
of the filament using these boundary densities. Accord¬ 
ingly, Ji, J p and J n are known for the entire filaments. 
For both (HL) 1 — (LL )2 phase region and (LL) 1 — (LL )2 
phase region, filament 2 is in LL phase and thus are 
known using Eq. ESI). Thus, the entire MF density and 
current profiles for all the species in both the filaments 
for the (HL) 1 — (LL) 2 phase region can be found out us¬ 
ing the known values of J\, J p and J n and P 2 L- For the 
(LL) 1 — (LL) 2 phase, pil is known using Eq. (1^51) . and 
along with J \, J p and J n are used to determine the den¬ 
sity and current profiles in this phase. The position of the 
shock of (+) particles in filament 1 (cc s i) is determined by 
matching the MF current solution of the (HL) 1 — (F/X) 2 
phase with (LL) 1 — (LL) 2 phase at the position of the 
shock, 

The position of the other shock on filament 2, x s2 , is 
determined by matching the MF current solution of the 
(HL) 1 — (HL) 2 phase with (HL) 1 — (LL) 2 phase at x s2 - 

J(HL) 1 -(LL) 2 ( Xs2 ') = J(HL)i — (HL) 2 (^s 2 ) (27) 

Fig.5 shows fairly good agreement for the density profiles 
derived from MC simulations and the MF predictions for 
the (HL) 1 — (HL) 2 and (LL) 1 — (LL) 2 - However for the 
(HL) 1 — (LT) 2 , the agreement between the MF solution 
and MC simulation does not match. 

This analysis has shown how inter-filament switching 
process leads to a wealth of new inhomogeneous phases, 
allowing for phase coexistence, and the possibility of 
shocks in both the filaments as opposed to the collective 
dynamics of transport in the absence of such filament 
interactions. 


B. Phase boundaries 


(SL) 1 — (SX) 2 : This arrangement, composed by three 
coexisting phases in the bulk, is characterized by shocks 
for (+)-particles in both filaments. In both filaments 
the right end is in (HL) 1 — (HL) 2 phase, while the re¬ 


in the previous subsection we have illustrated that any 
pair of pure phases are mediated by a phase coexistence 
region in the phase diagram. Whenever the current so¬ 
lutions corresponding to two different phases cross each 
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other at a particular spatial location along the filament, 
the system exhibits phase coexistence and the spatial lo¬ 
cation of the shock coincides with the location on the fil¬ 
ament where the two different current solutions intersect. 
The system selects the combination of those set of steady 
state density profiles for which the corresponding current 
is minimum at any given spatial location in the bulk. 
Exploiting this insight, we can now determine the entire 
phase diagram and the corresponding phase boundary by 
using the condition when the minimum value of the cur¬ 
rent along the filaments is allowed at one of their ends. 
Accordingly, the phase boundary separating any two re¬ 
gions in the phase space corresponds to the parameter 
values for which the location of the intersection of the 
different current solutions occurs at either of the filament 
ends. We obtain the MF phase boundary using the con¬ 
ditions for allowed phase in particular parameter range of 
entry and exit rate of particles and compare these results 
with MC simulations. We will concentrate on the phase 
diagram when varying {af — (3 +) while holding the other 
parameters constant. This choice clearly illustrates the 
qualitative new scenarios that filament switching brings 
into collective transport. Moreover, the scheme described 
can be straightforwardly extended to analyze the phase 
diagram when varying other sets of control parameters.. 

Boundary between {HL)\ — {LL) 2 phase with {HL) 1 — 
{SL) 2 '- Here the phase boundary is determined by the 
condition, 

J(hl) 1 -(ll) 2 ( 1 ) = ^{HL) 1 -(HL) 2 ^ ( 28 ) 

where, J(hl) 1 ~(ll) 2 (■*■) * s the current of (+) species in 
filament 2 at the right boundary at x = 1, when the sys¬ 
tem is in {HL) 1 — (LL ) 2 phase and J(h L ) 1 -(hl) 2 (■*■) 
the current of (+) species in filament 2 at x = 1 when 
the system is in {HL) 1 — {HL ) 2 phase. As discussed 
in the previous section, the {HL) 1 — {SL) 2 phase can 
be thought as an a mixture of the {HL) 1 — {LL )2 and 
{HL) 1 — {HL )2 phases. By matching the boundary con¬ 
ditions for {HL) 1 — {HL) 2 , we can determine J p , J n , 
Ji and at a; = 1, using Eq. (I2T1) for the 

boundary densities in {HL) 1 — {HL ) 2 phase. For the 
{HL) 1 — {LL )2 phase, we use Eq. (1751) in order to de¬ 
termine P 2 L > the boundary density for (+) at x = 0 in 
filament 2. Since there is overall particle conservation in 
the bulk, the values of J p , J n and J\ will be the same 
for both phases. Thus, using P 2 L along with the values 
of Jp, J n and Ji, we can now find out all the boundary 
values for densities in both filaments at x = 0 for the 
{HL) 1 — {LL )2 phase. Using the evolution equations for 
densities, Eqs. (116H19D . we can find out J‘^hl) 1 ~{ll) 2 at 
x = 1. Matching this expression for the current with 


J(hl) 1 -{hl) 2 f d x = 1 determines the MF phase bound¬ 
ary. 

Boundary between {HL)\ — {HL) 2 phase with {HL) 1 — 
{SL) 2 ~- Here the phase boundary is determined by the 
condition, 

J(HL) 1 -(LL) 2 (®') = ( 29 ) 

where, J^fjL)!-{LL)S 0) * s the current of (+) species in 
filament 2 at the left boundary when the system is in 
{HL) 1 — {LL) 2 phase and J(h L ) 1 -(hl) 2 (0) the current 
of (+) species in filament 2 at the left boundary when the 
system is in the {HL) 1 — {HL )2 phase. The procedure 
for finding out the phase boundary is same as in the 
previous case, except that the current matching is now 
done at x = 0. 

Boundary between {HL) 1 — {LH) 2 phase with {HL)\ — 
{LS) 2 '- The condition for phase coexistence in this case 
reads 

J(HL)i.-{LL)S^) = ( 30 ) 

where, J^hl) 1 -(ll) 2 (1) the current of (—) species in 
filament 2 at the right boundary when the system is in 
{HL) 1 — {LL) 2 phase and J(hl) 1 -(lh ) 2 (-0 i s the current 
of {—) species in lane—2 at the right boundary when the 
system is in {HL) 1 — {LH) 2 phase. For this case, J\ 
and J 2 can be determined since PiR,niR,p 2 L and n 2 L is 
known from Eq. (ED- Since for {HL) 1 — {LL )2 phase, n 2 R 
is known using Eq. m and n\R is already known using 
Eq. (1711) , thus J„ can be determined and hence using the 
known values of J\ and J 2 , J v can also be determined. 
This information is sufficient to determine the density 
and current profiles for both LL and LH phase in fila¬ 
ment 2. Therefore, the condition of current matching of 
(—) species on filament 2 at x = 1 determines the location 
of the phase boundary. 

Boundary between {HL) 1 — {LL) 2 phase with {HL) 1 — 
{LS) 2 '. The condition for this phase boundary for this 
case is, 

J (HL) 1 -(LL ) 2 (.°) = J (HL) 1 -(LH) 2 ^ ( 31 ) 

where, J(HL)!-(LL) 2 (ty ^he current of (—) species in 
filament 2 at the left boundary when the system is in 
{HL) 1 — {LL ) 2 phase and J(H L ) 1 -(LH) 2 (ty the current 
of (—) species in filament 2 at the left boundary when 
the system is in {HL) 1 — {LH ) 2 phase. The procedure to 
find the density and the current profiles follows exactly 
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the arguments as in the previous case. However, the 
condition of current matching of (—) species on filament 2 
is done at x = 0 and this determines the phase boundary 
between the two phases. 

Boundary between (LL )i — (LL) 2 phase with (SL)i — 
(LL) 2 - The condition to determine the phase boundary 
is given by 


^( 1 LL)i-(LL) 2 ( 1 ) — J(HL) 1 -(LL) 2 i 1 ) ( 32 ) 

Here, J(l L ) 1 -(ll) 2 ( 1) the curren t of (+) species in fil¬ 

ament 1 at the right boundary when the system is in 
(LL)\ — (LL) 2 phase and J(h L ) 1 -(ll) 2 ( 4 ) the current °f 
(+) species in filament 1 at the right boundary when the 
system is in (HL) 1 — (LL) 2 phase. For the (LL) 1 — (LL) 2 
phase the densities - p\R and p 2 R is known from Eq. 425D, 
which identifies J p . For the (HL) 1 — (LL) 2 phase pm, 
and n 2 R are determined by Eq. m and this is used 
to find Ji and J n . Having obtained these fluxes, the en¬ 
tire density and current profile for the two phases can 
be found out and the current matching condition for the 
(+) species in filament 1 determines the phase boundary 
in this case. 

Boundary between (HL)\ — (LL) 2 phase with (SL) 1 — 
(LL) 2 : Condition for this boundary in terms of boundary 
currents is given by 

^(LL)i-(LL) 2 (0) — T( 1 f ) : i ) 1 _( ii ) 2 ( 0) (33) 


is satisfied, where J(hr) 1 -(ll) 2 ^') the current of (+) 
species in filament 2 at the right boundary when the sys¬ 
tem is in the (HL) 1 — (LL) 2 phase and J^hr) 1 -(hr) 2 ( 1) 
is the current of (+) species in filament 2 at the right 
boundary when the system is in the (HL) 1 — (HL) 2 
phase. 

At the right boundary the system is in the (HL) 1 — (HL) 2 
phase, and therefore piR,P 2 R,nm and n\R can be de¬ 
termined using Eq. m- Consequently ,J p , J n and J\ 
are known and J^rl) 1 -(hl) 2 ^ can he found out. Since 
the left end of the filament 2 is in the LL phase, pir is 
known, and the entire density and current profile can be 
determined. From these profiles we identify the remain¬ 
ing flux, ■/^ i ) 1 _( i L) 2 ( 1)' Matching the current for the 
two profiles at the right boundary determines the phase 
boundary. 
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Here, J(r L ) 1 _( LL ) 2 ( 0) is the current of (+) species in 
filament 1 at the left boundary when the system is in 
(LL)i — (LL) 2 phase and J(j IL ) 1 -(ll)S 0) * s the curren f 
of (+) species in filament 1 at left boundary when the 
system is in (HL) 1 — (LL) 2 phase. For the (LL) 1 — (LL) 2 
phase the densities - Pir,Uir,p 2 r and n 2 R are known 
from Eq. (l25l) . which identifies J p and J n . For the 
(HL) 1 — (HL) 2 phase, p\R and n\R are determined from 
Eq. (l2lT) and this is used to find J\. Having obtained 
these fluxes, the entire density and current profile for the 
two phases can be found out and the current matching 
condition for the (+) species in filament 1 at the left 
boundary determines the phase boundary. 


Boundary between (SL) 1 — (LL)- 2 phase with (SL)i — 
(SL) 2 : The phase boundary derives from the parameters 
for which, 

'^(ffL)i-(LL) 2 ( 1 ) = ^ r (HL) 1 -(HL) 2 i l ) ( 34 ) 


FIG. 6: Phase space cut along a+—. Here, oq = 0.4, 
= 0.3, a+ = 0.8, &2 = 0.2, & = 0.25, 0 2 = 0.7 and 
7r = 1.0. (a) gives the MF phase diagram while (b) is obtained 
by MC simulation with N = 5000. 


Boundary between (HL)\ — (SL) 2 phase with (SL)i — 
(SL) 2 : The phase boundary is derived from the condition 

^(LL)i-(LL) 2 (°) = ^(ffL)i-(LL) 2 (°)’ ( 35 ) 

where is the current of (+) species in fil¬ 

ament 1 at the left boundary when the system is in the 
(LL)i — (LL) 2 phase and J(hr) 1 -(rr)S 1) * s t ^ ie current 
of (+) species in filament 2 at the left boundary when 
the system is in the (HL) 1 — (LL) 2 phase. 

At the right boundary the system is in the (HL) 1 — (HL) 2 
phase. Therefore pir, p 2 R , U\r and n 2 R can be de¬ 
termined using Eq. (ED. while n 2 R can be determined 
using Eq. (l25l) . which identifies J p , J n and J\. Using 
this information the entire density and current profile 
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FIG. 7: Phase space cut along a+— dj 1 " plane. Here, a 1 = 0.2, 
df = 0.7, aj = 0.2, = 0.2, d^ = 0.2, /3J = 0.6 and 

7r = 1.0. (a) gives the MF phase diagram while (6) is obtained 
by MC simulation with N = 5000. 


can be determined for the phases (LL )i — (LL) 2 and 
(HL)\ — (LL) 2 - Matching the current for the two pro¬ 
files in these phases at the left boundary determines the 
phase boundary. 


LL phase to and SL phase with a shock developing on 
this other filament. The phase diagram also shows the 
possibility to have shock reentrant phases. For a given 
entry rate, the increase in the exit rate naturally favours 
a transition from HL to LL phases, but these transitions 
are modulated by the developments of shocks. As a re¬ 
sult, in the transition from HL-LL to LL-LL is mediated 
by an intermediate region of SL phases, and for large 
enough entry rates, the LL phase is destabilized by the 
development of an SL phase as the exit rate increases. 



All the other possible combination of phases bound¬ 
aries can be obtained by simply interchanging the labels 
of filament 1 and 2 and using exactly the same condi¬ 
tions for phase boundaries that have been described in 
this subsection. 

From all these conditions, we can now build a complete 
phase diagram. Fig.6 and Fig.7 show the comparison of 
MF and MC phase diagram in different phase plane cuts, 
as a function of (a] 1- — /?i~) for fixed values of the rest of 
the control parameters. The MF phase boundaries that 
have been obtained exhibit fairly good agreement with 
the phase diagrams obtained by MC simulations, which 
shows that MF is quantitatively accurate to describe the 
different phases that characterize transport intros sys¬ 
tem. The contrast between Fig.6 and Fig.7 also shows 
that the topology of the phase plane can drastically be 
altered by tuning the parameters corresponding to the 
particle entry and exit rates although the filaments them¬ 
selves are weakly coupled through the filament switching 
processes. Further, we also see that changing the particle 
entry and exit rates in one filament can affect the phases 
in the neighbouring filament. For instance in Fig.7 as one 
increases the value of entry rate of (+) particles in fila¬ 
ment 1, keeping the exit rate of (+) particles fixed, the 
resultant phase in the other filament passes over from an 


FIG. 8: Phase plane cut along a~t~ di" plane, [(a) and (c)] are 
without filament switching (n = 0) as discussed in 0,a. m 
and (d)] are with filament switching (n = 1). In (&); aq = 0.4, 
/If = 0.3, at = 0.8, aq = 0.2, fit = °- 25 , Ps = 0.7. In (d); 
or = 0.2, di“ = 0.7, at = 0.2, aq = 0.2, fit = 0.2, dr = 0.6. 
Parameter regime of (a) is same as that of (c) and parameter 
regime of (b) is same as that of (d). 


The resulting phase diagram for this system can topo¬ 
logically be very distinct from the phase diagrams for 
two species transport in the absence of particle switch¬ 
ing between filaments. To emphasize this fact, in Fig. 8 
shows the MF phase diagram for a system in the absence 
of filament switching (n = 0), Figs.8a and 8c, with cor¬ 
responding predictions when particle filament exchange 
is allowed Figs.8b and 8d keep‘the same weak exchange 
rates between the two filaments and modify the corre¬ 
spond to a choice of entry and exit rates. We have chosen 
representative parameters to show the differences associ¬ 
ated to particle filament exchange. For Figs.8a and 8b 
one illustrates that one of the new features introduced 
by particle exchange between filaments is the appear¬ 
ance of phase coexistence regions sandwiched between 
pure phases; a scenario forbidden when opposite parti¬ 
cles displace along a unique filament [T, 2|. Figs. 8(c) 
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and Fig.8(d), for a different set of parameter values, show 
that particle filament exchange can have a deeper im¬ 
pact on the phase diagram topology. In this case we see 
that, while in the absence of filament switching dynamics 
the system always in a pure phase, the effect of filament 
switching of particles manifests as enriching the phase 
behaviour for the system which allows for phase coexis¬ 
tence and presence of shocks in the two filaments. 

As illustrated in these two examples, generically we 
find that two pure phases are always connected by a 
phase region where the system exhibits phase coexistence 
and bulk localized shocks. 


V. CONCLUSIONS 


To summarize, we have studied a multi-filament driven 
system with oppositely directed species of particles when 
the filaments are weakly coupled. Particle filament 
switching processes constitute correlated events because 
particles can only swap filaments, with a prescribed fi¬ 
nite probability, when oppositely directed particles meet 
on the same filament. This aspect of filament switching 
mimics cellular cargo switching between neighbouring fil¬ 
aments during intracellular transport. We find that the 
interplay of the entry and exit processes of particles at 
the filament boundaries has a profound impact in the 
collective organization of the two species of displacing 
particles, leading to a variety of new scenarios. Specifi¬ 
cally, we have identified the development of phase coex¬ 
istence on the filaments, inhomogeneous density profiles, 
density shocks localized in the bulk and bidirectional cur¬ 
rent flows in the system. We have developed a mean field 
theory (MF) to characterize these phenomena, and have 
shown that the steady state density and current profiles 
of particles and the phase diagram obtained using a MF 
formulation match reasonably well with the Monte Carlo 
(MC) simulation results. 

While in this paper we have focused on the implica¬ 
tions that weak coupling between filaments would have 
on transport when there are particle input and output 
in both filaments, it would be interesting to explore the 
regime of strong filament coupling, where we expect a 
weaker spatial inhomogeneity in the particle density pro¬ 
files. Further, for many biological situations such as 
transport in axons, it is a priori not clear whether bound¬ 
ary loading and off-loading of cellular cargo happens for 
all the parallel filaments or for specific filaments. In such 


situation one needs to determine the steady state distri¬ 
bution of cargoes and the resultant phases when one of 
the filament has much higher particle entry rates than 
the other. 

As an extreme case, we have considered a situation 
where one filament has closed boundaries. Starting from 
a random distribution of particles in both filaments, we 
have observed the development of phase segregation be¬ 
tween (+) and (—) particles in the closed filament. This 
phenomenon happens only due to the correlated lane 
switching process. The resulting phase segregated state 
in the blocked filament does not have any flux. Fol¬ 
lowing the time evolution of such a system shows that 
starting from a random configuration, all the vacancies 
are expelled and eventually the filament stops exchanging 
particles with the other filament as the blocked filament 
attains a jammed configuration, with the (+) particles 
piling up from the left end of the filament while the (—) 
particles pile up on the right end . Understanding the 
transition of this phase segregated jammed steady state 
to the steady states discussed in this paper as one slowly 
increases particle input and output for the filament which 
is initially closed at the boundaries, remains as an inter¬ 
esting open challenge. 


Appendix A: Numerical simulations 


For determining the steady state density and current 
profiles on filaments, Monte Carlo (MC) simulations have 
been performed to simulate the various processes de¬ 
scribed in Section[Tl] For a MC move, a filament is chosen 
at random and then a site in that particular filament is 
chosen at random with equal probability. If a particle 
(+) or (—) is present then a move is made for the various 
dynamic processes (e.g; translation or lane switching), 
proportional to the respective rates. We begin the simu¬ 
lation run starting from a random initial distribution of 
particles in the two filaments and let the system evolve 
and reach steady state. After the system has attained 
steady state, averaging is done over occupation number 
and current in the lattice. Typically we wait for initial 
transient of 1000^ swaps, where q is rate of the slowest 
process among all the different dynamic processes occur¬ 
ring in the lattices. We have further checked that the 
system indeed reaches its steady state by comparing the 
final density and current profiles at the end of the tran¬ 
sient. We then collect the data for occupation number 
and current with a period > 10y and average them over 
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5000 time swaps. 

In order to determine the phase boundaries by MC 
simulations, we use the fact that all phase transitions 
between pure phases are mediated by phase coexistence 
regions with shocks in the density profile. The phase 
boundaries can be determined numerically by tracking 
when the location of shock reaches the filament boundary. 
However, due to the finite size effects of the system, the 
shock has a certain finite width. We determine the shock 
width and track the position of the midpoint of the shock 
to identify the phase boundaries numerically and decide 
when they have reached a filament end. We have used 
system size of N = 5000 for determining the location 
of phase boundaries. For a fixed /3+ we have varied 
in steps of 10 -3 and this sets the accuracy of the phase 
boundaries obtained numerically. 


Appendix B: Choice of branches for the MF solution 
for density 


To choose a proper branch uniquely from the poten¬ 
tial eight solutions that can be derived from the Mean 
Field solutions for a particular variable we have to look 
carefully at the density profiles of various species in a 
particular phase. In Eas. (I161UT)1) . notice that various 
combinations of ry^’s, yt^’s and z> ±! s appear in this set 
of differential equations which govern the spatial den¬ 
sity profile for each of the individual species. Each of 
the individual combinations of ry^’s, yt^’s and ap¬ 
pearing in these set of differential equations are exclu¬ 
sively functions of a single variable. For example if we 
consider ryj^, then it appears in Eq. (flbl) , where the ex¬ 
plicit form reads as ?y+ = | + \J\ +f>i(l — Pi) — Ji and 

7 bi = \ - \]\ +Pi ( 1 -Pi) - Ji- Each jy+ and ?y“ could 
separately be combined with each of the two different 
values of ytj^ and that appear in the expressions of 
Eq. (IBID . As there are eight possible combinations of ry±, 
yt^ and , therefore there would be eight potential solu¬ 
tions corresponding to the choice of a particular bound¬ 
ary condition. Similarly there would be eight possible 
solutions for Eas. (11711191) . 

To do the classification we have to look at the various 


expressions for 77=*= ’s, yt ±: s and z^’s and they read as, 

Vp! = \ ± \j\ +Pi(l -Pi) - Ji 
hp, = \ ± \j\ +Pi(l ~Pi) ~ J P 
"p! = \ ± \j{\ ~ Pi) 2 + Ci (Bl) 


Vnt = \±^\+ni(l-m)-Ji 
= \ ± y \ +ni{l -ni) - J„ 

= \ ± \j{\ -m) 2 -C 2 (B2) 


Vp 2 = \ ± \j\+P2(l - P2) - Ji 
Pp 2 = \ ± \]\+P^-P 2 )-Jp 
4 = \ ±^{\~P2) 2 + C 2 (B3) 


vt 2 = 7y± \J\+n 2 {\ - n 2 ) - J 2 

Pn 2 = \ ± \j\+n 2 {l-n 2 ) ~ J n 

v n 2 = \ ± “Ci (B4) 

Among the eight different branches of Ea. fTGl) which 
arise due to eight possible combinations of ry±, yt^ and 
v ± x , we take that particular branch of the equation which 
is a combination of 7 )~ x , ii~ x and . This branch would 
be classified as Sol—1. The entire nomenclature is clas¬ 
sified in Table-U for different solutions of Eq. (flbl) . The 
same nomenclature is true for the different solutions of 
Eqs. (fT7!([H)|) . 

Let us assume that filament 1 is in LL phase whereas 
filament 2 is in HL phase. Therefore the density of (+) 
and (—) particles in both the filaments and hence the 
variables will have the following values. 

1 1 

Pl< 2’ ni< 2 

P2 > \, n2< \ 
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TABLE I: Classification of the branches 


V 

V 

V 

Branch 

Vri 

Vpi 

v p-\ 

Sol-1 

Vp-i 

Vp-i 

JS+ 

u p 1 

Sol-2 

V, 

1 1 P1 

V P1 

Sol-3 

Vp, 

Vpi 

V P1 

Sol-4 

Vp! 

V'pi 

U V1 

Sol-5 

Vpi 

Vpi 


Sol-6 

Vp, 

V'pi 

U P1 

Sol-7 

Vpi 

Vpi 

"pi 

Sol-8 


Thus, the proper choice of roots, in this case, while sub¬ 
stituting other variables in terms of pi would be 


m = Ji = Vp, 

Pi = \ + \j\ +Pi(l -Pi) - J P = Mpi 
m = \ - \j{\ -Pi) 2 3 4 5 6 * 8 9 * 11 12 13 14 + Ci = v~ (B5) 

Thus sol—3 is the correct branch for pi when the system 


is in (LL)i — {HL )2 phase. Following the same procedure 
we can show that in this particular phase sol—2, sol—1 
and sol—4 are the proper choices for the variables n 1 , 
P 2 and n 2 respectively. Relevant branches for all other 
phases are depicted in Table-HIl 


TABLE II: MF Solutions In Different Phases 


Phase 

Pi 

m 

P2 

n 2 

(LL) 1 - {LL) 2 

Sol-1 

Sol-1 

Sol-1 

Sol-1 

{LL) 1 - (HL) 2 

Sol-3 

Sol-2 

Sol-1 

Sol-4 

{LL)i - ( LH) 2 

Sol-2 

Sol-3 

Sol-4 

Sol-1 

{HL) 1 - {HL) 2 

Sol-3 

Sol-6 

Sol-3 

Sol-6 

{HL) 1 - (LH) a 

Sol-2 

Sol-7 

Sol-7 

Sol-2 

{LH) 1 - {LH) 2 

Sol-6 

Sol-3 

Sol-6 

Sol-3 
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